import numpy as np, vtk , cv , thread

imsize = (256,256)

def ReCastFloatToUShort(data):
    print np.can_cast(np.float32, np.uint16)


#Costas 6/29/2011
def LThresholdData(thresh,data):
    data[data<thresh]=0
    return data

def UThresholdData(thresh,data):
    data[data>thresh]=0
    return data



def Reconstruct(data):

    splatter = vtk.vtkGaussianSplatter();
    splatter.SetInput(data.GetOutput());
    splatter.SetSampleDimensions(50,50,50);
    splatter.SetRadius(0.5);
    splatter.ScalarWarpingOff();
    return splatter
 
def OrthoRayCast(data,thres):
    a=0.0
    dataf=data
    
    img=cv.CreateImage(imsize,cv.IPL_DEPTH_8U, 1);
    cv.SetZero(img)
    print data.shape
    found=False
    for i in range(data.shape[0]):
        for j in range(data.shape[2]):
            found=False
            foundB=False
            for k in range(data.shape[1]):
                if(dataf[i,k,j]>thres and not found):
                    cv.Set2D(img,imsize[1]-data.shape[2] + j,imsize[0]-data.shape[0] + i,data.shape[1]-k-1);
                    #dataf[i,k,j]=k
                    found=True
                    break;		

        cv.ShowImage('depthmap',img)
        cv.SaveImage('b.png',img)	

    return dataf

def ConvertContourToArray(contour):
    pointSampler = vtkPolyDataPointSampler();
    pointSampler.SetDistance(.1);
    pointSampler.SetInputConnection(contour.GetOutputPort());
    pointSampler.Update();
    return pointSampler()

def FullOrthoCast(data,thres,save,ThinShell):

     
    #fp=open('model.csv','w')
    a=0.0
    dataf=np.zeros(data.shape, dtype=np.uint8)
    print "executing FullOrthoCast"
    
    img=cv.CreateImage((data.shape[0],data.shape[2]),cv.IPL_DEPTH_8U, 1);
    cv.SetZero(img)
    for i in range(data.shape[0]/4,3*data.shape[0]/4):
        for j in range(data.shape[2]):
            found=False
            foundB=False
            for k in range(data.shape[1]/2):
                if(data[i,k,j]>thres and not found):
                    if(save):
                        pass#fp.write(str(j)+','+str(i)+','+str(k)+'\n')					
                    cv.Set2D(img,j,i,k);
                    dataf[i,k,j]=300
                    found=True

                elif found and ThinShell:
                    dataf[i,k,j]=300
            for k in range(data.shape[1]-1,data.shape[2]/2-10,-1):
                if(data[i,k,j]>thres and not foundB):
                    if(save):
                        pass#fp.write(str(j)+','+str(i)+','+str(k)+'\n')		
                    dataf[i,k,j]=300
                    foundB=True
                elif  foundB and ThinShell:
                    dataf[i,k,j]=300

    for j in range(data.shape[2]):
        for k in range(data.shape[1]):
            foundR=False
            foundL=False
            for i in range(0,data.shape[0]/4):
                if(data[i,k,j]>thres and not foundR):
                    if(save):
                        pass#fp.write(str(j)+','+str(i)+','+str(k)+'\n')		
                    dataf[i,k,j]=300
                    foundR=True
                elif  foundR and ThinShell:
                    dataf[i,k,j]=300
            for i in range(data.shape[0]-1,3*data.shape[0]/4-10,-1):
                if(data[i,k,j]>thres and not foundL):
                    if(save):
                        pass#fp.write(str(j)+','+str(i)+','+str(k)+'\n')		
                    dataf[i,k,j]=300
                    foundL=True
                elif foundL and ThinShell:
                    dataf[i,k,j]=300
    #fp.close()
    return dataf
def FrontOrthoRayCastNew(data,thres,name):
    """Compute coronal depath map
    data: 3D numpy array representing original volume
    thres: threshold for which values less than thres are ignored in generating depth map
    name: filename for file
    """
    mask = np.where(data>thres,1,0)
    depthMap = mask*(data.shape[1]-np.argmax(mask,axis=1)).astype(np.uint8)
    kmax = depthMap.max()
    tmp = os.path.splitext(name)
    img = cv.GetImage(cv.fromarray(depthMap))
    minV,maxV,min_loc,max_loc = cv.MinMaxLoc(img)
    cv.SaveImage(tmp[0]+'_coronal.png',img)		
    return img,maxV,max_loc

def FrontOrthoRayCast(data,thres,name):
    """Compute coronal depath map
    data: 3D numpy array representing original volume
    thres: threshold for which values less than thres are ignored in generating depth map
    name: filename for file
    """
    a=0.0
    dataf=np.zeros(data.shape, dtype=np.uint8)
    print "executing DepthImageGeneration"
    
    img=cv.CreateImage(imsize,cv.IPL_DEPTH_8U, 1);
    cv.SetZero(img)
    print data.shape
    found=False
    kmax  = 0;
    max_loc =(0,0)
    for i in range(data.shape[0]):
        for j in range(data.shape[2]):
            found=False
            foundB=False
            for k in range(data.shape[1]/2):
                if(data[i,k,j]>thres and not found):
                    cv.Set2D(img,(imsize[1]-data.shape[2])/2 + j,(imsize[0]-data.shape[0])/2 + i,data.shape[1]-k-1);
                    if(kmax<data.shape[1]-k-1):
                        kmax = data.shape[1]-k-1
                        max_loc = ((imsize[0]-data.shape[0])/2 + i,(imsize[1]-data.shape[2])/2 + j)
                    break;

    print kmax				
    for i in range(data.shape[0]):
        for j in range(data.shape[2]):
            g=cv.Get2D(img,(imsize[1]-data.shape[2])/2 + j,(imsize[0]-data.shape[0])/2 + i)
            if(g[0]>0):
                cv.Set2D(img,(imsize[1]-data.shape[2])/2 + j,(imsize[0]-data.shape[0])/2 + i,g[0]-(data.shape[1]-kmax));


    # cv.Smooth(img, img, cv.CV_MEDIAN, 3)	
    print "finished DepthImageGeneration"
    # GoodFeatures(img)				
    #cv.ShowImage('depthmap',img)
    
    saveloc= name[:len(name)-4]
    
    cv.SaveImage(saveloc+'.png',img)		
    # cv.WaitKey(0)		
    return img,kmax,max_loc

def GoodFeatures(img):
    print "started Points"
    im2 = cv.CloneImage(img)
    eigimg=cv.CloneImage(img)
    tempimg=cv.CloneImage(img)
    points = cv.GoodFeaturesToTrack(img,eigimg,tempimg,100, 0.1,.2)
    print points
    for p in points:
        cv.Circle(im2, (int(p[0]),int(p[1])) , 5,(0,255,0), 1)
    
    cv.ShowImage('points',im2)

def ClipSurface(surface):
    plane=vtk.vtkPlane()
    plane.SetOrigin(0,0,0)
    plane.SetNormal(1,1,0)
    
    extract = vtk.vtkExtractPolyDataGeometry();
    
    extract.SetInput(surface.GetOutput());
    
    extract.SetImplicitFunction(plane);
    
    extract.ExtractInsideOff();
    
    extract.ExtractBoundaryCellsOn();
    
    extract.Update();
    return extract
   
def ReconstructSurface(img):
    surf = vtk.vtkSurfaceReconstructionFilter()
    surf.SetInputConnection(img.GetOutputPort())
    return surf

def ExtractSurface(volume):
    #Costas 6/29/2010
    contours = vtk.vtkContourFilter()
    contours.SetInput(volume.GetOutput())
    #contours.SetValue(0,10000)
    contours.GenerateValues(1,10,11)
    contours.SetUseScalarTree(1)
    # contours.SetNumberOfContours(10)
    return contours

def Decimate(surface,level):
    Deci =vtk.vtkDecimatePro()
    Deci.SetInput(surface.GetOutput())
    Deci.SetTargetReduction(level)
    
    #deci.PreserveTopologyOn()
    return Deci


def ExtractLargest(volume):
    confilter =vtk.vtkPolyDataConnectivityFilter()
    confilter.SetInputConnection(volume.GetOutputPort())
    #confilter.
    confilter.SetExtractionModeToLargestRegion()
    
    return confilter

def TriangleStrip(volume):
    Stripper = vtk.vtkStripper()
    Stripper.SetInputConnection(volume.GetOutputPort())
    normals = vtk.vtkPolyDataNormals()
    normals.SetInput(Stripper.GetOutput())
    normals.FlipNormalsOn()
    
    return normals

def SmoothSurface(surface,noIter):
    smoother = vtk.vtkSmoothPolyDataFilter()
    smoother.SetInput(surface.GetOutput())
    smoother.SetNumberOfIterations(noIter)
    normals = vtk.vtkPolyDataNormals()
    normals.SetInput(smoother.GetOutput())
    normals.FlipNormalsOn()
    
    return normals




def GetOuterSurface(data,renderer):

    poly = vtk.vtkImageDataGeometryFilter()
    poly.SetInputConnection(data.GetOutputPort())
    poly.Update()
    
    
    polydata = poly.GetOutput()
    
    silhouette =  vtk.vtkPolyDataSilhouette()
    
    silhouette.SetInput(polyData);
    silhouette.SetCamera(renderer.GetActiveCamera());
    silhouette.SetEnableFeatureAngle(0);
    
    return silhouette
    

def GetVOI(volume):
    extract = vtk.vtkExtractGrid()
    extract.SetVOI((0,10,0,10,0,10))
    extract.SetInputConnection(volume.GetOutputPort())
    return extract



def MarchingCubes(volume):
    surface = vtk.vtkMarchingCubes();
    surface.SetInput(volume.GetOutput())
    surface.ComputeNormalsOn()
    surface.SetValue(0, 1000);
    return surface
    
    

